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Abstract 

Supernova remnants (SNRs) are one of the most energetic astrophysical events and are 
thought to be the dominant source of Galactic cosmic rays (CRs). A recent report 
on observations from the Fermi satellite has shown a signature of pion decay in the 
gamma-ray spectra of SNRs. This provides strong evidence that high-energy protons are 
accelerated in SNRs. The actual gamma-ray emission from pion decay should depend 
on the diffusion of CRs in the interstellar medium. In order to quantitatively analyse 
the diffusion of high-energy CRs from acceleration sites, we have performed test particle 
numerical simulations of CR protons using a three-dimensional magnetohydrodynamics 
(MHD) simulation of an interstellar medium swept-up by a blast wave. We analyse the 
diffusion of CRs at a length scale of order a few pc in our simulated SNR, and find the 
diffusion of CRs is precisely described by a Bohm diffusion, which is required for efficient 
acceleration at least for particles with energies above 30 TeV for a realistic interstellar 
medium. Although we find the possibility of a superdiffusive process (travel distance 
oc in our simulations, its effect on CR diffusion at the length scale of the turbulence 
in the SNR is limited. 
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1. Introduction 

Supernova remnants (SNRs) have long been believed to be the source of hadronic 
Galactic cosmic rays (GCRs) up to energies of the ‘knee’, near 5 x 10^® eV, of the 
cosmic ray (CR) spectrum. Supernova explosions forming collisionless shock waves induce 
the shocked gas and relativistic particles (hereafter cosmic rays) that produce multi¬ 
wavelength thermal and nonthermal emission. Diffusive shock acceleration (DSA) is the 
most promising mechanism for converting the kinetic energy of a supernova explosion 
into energetic particles Biii and plays an important role in nonthermal emission 
during the overall process, (e.g., Q). 
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In the framework of DSA, an individual charged particle experiences many collisions 
with background electromagnetic waves and gains energy by shock crossing. This leads 
to a nonthermal CR spectrum of the power-law form N{e) ^ e“^. Shock acceleration 
by DSA in SNR shocks is associated with transport processes, and some of the highest 
energy CRs eventually escape from their acceleration sites by a so-called diffusion process 
due to interactions with turbulent magnetic fields. To analyse the diffusion process we 
have to determine the effective diffusion coefficient D{r,p), and the understanding of 
D(r,p) is necessary to interpret many astronomical observations. 

Several studies of escaping CRs have already been done and a strong spatial corre¬ 
lation between TeV emission and the molecular gas distribution at the Galactic Center 
has been observed Sii- The pion-decay signature in SNRs is believed to be evi¬ 
dence for protons accelerated in middle-aged SNRs interacting with molecular clouds, 


[E El, ESQ EH- Recent observations suggest that gamma-ray emission and 
CO-1-HI emission are spatially correlated in young SNRs RX J1713.7-3946 and RX 
JOS52.0-4622 (EE , see also EE El. 

Nevertheless, the identification of pion-decay gamma rays is difficult because high- 
energy electrons also produce gamma rays via bremsstrahlung and Inverse Compton 
scattering (leptonic model) El- X-ray observations show that electrons are accelerated 
to highly relativistic energies in SNR shocks El- Therefore, in order to understand 
the acceleration sites of CRs, it is crucial to distinguish GeV-TeV emission from Inverse 
Gompton scattering by CR electrons and the decay of neutral pions produced by inelastic 
collisions between CR protons and ambient thermal nuclei. 

In this paper, we investigate the diffusion of CRs using a hydrodynamics simulation 
of a strong shock wave propagating in a realistic multiphase interstellar medium and a 
one-phase medium. The organization of the paper is as follows. In Section 2, we describe 
the three-dimensional hydrodynamics simulations and the resulting configuration of elec¬ 
tromagnetic field used. We also briefly introduce the process of Bohm Diffusion. The 
results of test particle simulations performed in these environments are shown in Section 
3. We investigate the properties of escaping CRs in terms of the diffusion coefficient in 
both energy and configuration space, and finally we summarize and discuss our findings 
in Section 4. 


2. Methods 

2.1. Setup of Background Interstellar Medium 

2.1.1. Multiphase medium 

The interstellar medium (ISM) is an open system in which radiative cooling and 
heating are effective. It is an inhomogeneous, multiphase system in which gases of differ¬ 
ent temperatures, densities, and ionization fractions can coexist in approximate pressure 
equilibrium. Diffuse warm gas (diffuse intercloud gas) with T ~ 10"* K and HI clouds 
(interstellar clouds) with T ~ 10^ K are approximately in pressure equilibrium in a typ¬ 
ical ISM environment. As a consequence of the thermal instability driven by external 
compressional events such as shock waves due to expanding HII regions or very late 
phase SNRs, unstable gas evolves into diffuse gas and HI clouds E, El, E- Therefore, 
inhomogeneities inevitably emerge and remain ubiquitous in the ISM. The characteristic 
length scale of an inhomogeneity can be expressed in terms of the “Field length", which 
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is the critical wavelength of the thermal instability [33|, [3^. The Field length depends 
on density and temperature and can be smaller than 1 pc. A blast wave generated by 
supernova expansion sweeps up the dense and clumpy HI clouds of the multiphase ISM, 
which eventually generates strong velocity shear in the magnetic fields. Magnetic fields 
undergo amplification from thei r typica l strength of /rG to mG due to the turbulent dy¬ 
namo in the post-shock region ( [35|, 127], [33 [STf ). This process may explain the existence 
of magnetic helds of mG strength investigated by Y. Uchiyama et al. (2007) [s^. 

T. Inoue et al. (2012) performed ideal three-dimensional magnetohydrodynamic 
(MHD) simulations of a strong shock wave {vsh ~ 2500 km s“^) propagating in a re¬ 
alistic multiphase ISM as the pre-shock region. We use the data for the perpendicular 
shock of T. Inoue et al. (2012) at t = 750 years as the background ISM to set up the 
electromagnetic field for our microscopic particle simulations. 

To generate a multiphase ISM, T. Inoue et al. (2012) solved the ideal MHD equations 
including cooling, heating, and thermal conduction, which determine the unstable scale of 
thermal instability. They considered a net cooling function and photoelectric heating, and 
generated an inhomogeneous medium via thermal instability. The simulation considered 
ideal gas and used an adiabatic index of F = 5/3. The mean number density, initial 
thermal pressure, and initial magnetic field strength were taken to be (no) = 2.0 cm“^, 
p/ks = 2887 K cm“^, and Boy = 5.0 pG, respectively, at the x = 0 boundary plane. For 
the density, they imposed random density fluctuations with a thermally unstable state 
in the range 10 K < T < 10'^ K for effective cooling and heating. In the resulting clumpy 
cloud, they induced a high Mach number shock wave using a hot plasma with Ph/ks — 
10® K cm“^ and (uh) = 0.1 cm“^. 


2.1.2. One-phase medium 

T. Inoue et al. (2013) performed ideal three-dimensional magnetohydrodynamic 
(MHD) simulations to investigate the interaction between blast wave (vsh 1800 km s“^) 
and interstellar density fluctuations. They investigated the magnetic field amplification 
and the magnetic field distribution of turbulent SNRs driven by the Richtmyer-Meshkov 
instability (RMI). They assumed an adiabatic gas with adiabatic index r=5/3 and used a 
high Mach number shock wave. Density fluctuations superposed by sinusoidal functions 
were included and followed an isotropic power-law spectrum with random phases. The 
power spectrum of the density fluctuations was shown to be described by an isotropic 
power law for the wavenumber k in the inertial range of turbulence: P{k) = p\k^ oc , 
where pk is the Fourier component of the density. The mean number density, the initial 
thermal pressure, and the initial magnetic held strength were taken to be (no) = 0.5 
cm“^, p/ks =4 X 10^ K cm“^, and Bn = 3.0 p,G, respectively. The parameters repre¬ 


sent typical values in the diffuse ISM 4^ . To induce the blast wave, a hot plasma is 
set up as follows: ph/ks = 2 x 10® K cm“®, (n^) = 0.05 cm“®, and Boy = 3.0 pG at 
the a: = 0 boundary plane. This creates the primary shock wave whose normal vector 
is perpendicular to the mean magnetic held. The SNR is modelled as a young SNR 
(age = 10® years) with a late with velocity of 1800 km s“^. We used the data for the 
perpendicular shock on model 1 in at t = 700 years as the background ISM. 

A simulation box with size L^ox = 2 pc is used, and the system resolution is Ax = 
T6oa:/(number of grid cells) = 1.95 x 10“® pc, where the number of uniform grid cells 
is 1024®. Periodic boundary condition is used for the yz-plane. The above simulation 
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resolution and boundary condition are applied for both media. The Fourier pow er spect ra 
of SNR turbulence for both media are given in Figure [TJ (for more details see 


27, 


2.2. Test particle simulations 

The magnetized turbulent medium flows along the (positive) x-direction. We calcu¬ 
late the trajectories of CR particles using a snapshot of the MHD simulation data. We 
consider a collisionless environment for CRs because the mean free path, A = l/ntr, of 
relativistic particles is much larger than the gyroradius of CR particles. In the framework 
of non-relativistic ideal MHD, electromagnetic fields can be expressed as 


E{t, r) 


X B(t,r) 
c ’ 


( 1 ) 


where V is the fluid velocity, B is the local magnetic field at the position of the particles, 
and c is the speed of light. When we compute the Lorentz force acting on each CR 
particle, we interpolate V and B at the position of the particle before we calculate E 
using Equation O- To trace the trajectories of CR, we solve the relativistic equation 
of motion using a fourth-order Runge-Kutta method for updating particle information. 
The momentum of a particle is defined as 


p = jmiV, 


( 2 ) 


where v is the velocity of a particle. 7 is the particle Lorentz factor, and rrii is the mass 
of particle i. The equations governing the motion of relativistic charged particles are 
expressed as 


dp 

dt 

dr 

dt 


Q^{E + 


vx B 


P 

7771 *’ 


( 3 ) 

( 4 ) 


where qi is the charge of particle i. The gyroradius, Rg, of particle i can be expressed by 


R„ 


'yniiC" 


for 


( 5 ) 


The time discretization. At, for each particle is as follows: 


tr — Qg /N^ 


( 6 ) 


Ax 

■ max{\va:\,\vy\, |uz|) 


( 7 ) 


At = min(tr, tg), 


( 8 ) 


where fig = \ qi\B /jmiC denotes the particle gyrofrequency and Ax is the size of a grid cell. 
Nd is the number of divisions which determines the resolution of a particle trajectory. The 
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time discretization is determined by inverse of gyrofrequency as it shows in Equation 
To trace the particle trajectory, At should be smaller than gyrofrequency. is selected in 
our simulation when the particle energy is sufficiently small such that Rg < Ax, whereas 
ts is selected for higher energy cases satisfying Rg > Ax. For more details on Rg, see 
Figure [3] in Section 3.2. We obtain At by considering itg and the velocity components 
(vx, Vy, and Vz) of a particle as well as using a sufficiently large value of Nd for tracing 
the trajectories of particles. For example, At was about 250 s for eo = 10^ TeV. In all 
calculations, we used Nd = 10^, which is the value for keeping the errors of Runge-Kutta 
time-integration sufficiently small (< 10“^® eV) for energy conservation. 


3. Results of Numerical Simulations for Diffusion of Cosmic Rays 

To calculate the transport of particles, we released sets of 10^ protons with the same 
initial energies but isotropic velocity distributions in the computational frame. The range 
of the initial energies (eo) considered is 10^^ eV to 10^® eV at intervals of log Eq = 0.5 eV. 
We assigned particle positions randomly in a cubic box region in the post-shock region 
(see Tables 1 and 2). Initially the directions of velocities are also random. To obtain 
the background fluid velocity and magnetic field we applied a linear interpolation to 
the MHD simulation data because the magnetic and fluid fields are defined at discrete 
points. We use periodic boundary conditions in the j/z-direction as was done in the MHD 
simulation. We do not follow the time integration for orbits of particles that have escaped 
the a;-boundary. First, we show the results of a multiphase medium case. 


3.1. Diffusion process in a Multiphase Medium 

3.1.1. Diffusion in Energy Space 

In a turbulent medium, second-order Fermi acceleration is expected where some parti¬ 
cles gain energy stochastically for a sufficiently long time interval. This can be considered 
to be a random walk in energy space (Ae oc '/i). 

We calculated up to 10^ years for 10^ TeV and 10^ TeV protons. Figure [3] show the 
magnitude of the particle energy, Ae, as a function of logt. Dashed lines in both figures 
represent energy gain proportional to t^'^. This demonstrates the stochastic behaviour 
of energy gain in our simulation. However, the energy gain rate is very small within 10^ 
years, less than 1% of the initial energy. Thus, we describe the energy in this paper with 
initial energy (^o)- 


3.1.2. Diffusion in Configuration Space 

Gyroradius (Equation Q) in units of pc for different values of £o shows in Figure 
[31 The values are calculated for 10^ particles at five different values of Eq at intervals 
of logEo = 0.5 eV. Each symbol of the solid lines corresponds to the gyroradius at a 
different time. Black open and filled circles correspond to t = 10 years and 10^ years, 
respectively. Red stars indicate for an initial time. Note that the smallest (Tmm) and 
largest {Lmax) turbulent scales in our simulation correspond to a few grids and a few 
hundreds grids, respectively 23. The grid size corresponds to Ax = 6.03 x 10^® cm. 
Particles are supposed to be reflected by the turbulent field when Lmm < Rg < Lmax- 
For the case of Rg < Lmin-, particles only see an almost uniform field due to grid based 
method, and thus are expected to undergo E x B drift. In this case the travel distance 
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is proportional to t at early times. In Figure O particles with Eq = eV, 10^"^ eV, 

10i4-5 gY^ 2.0^® eV satisfy Lmin < Rg < Lmax- We do not consider the case of Eq = 
10^® eV for our study of diffusion process in configuration space, because the number of 
escaped particles in our simulation is large. 

To analyze the diffusion in configuration space we measure the displacement of trans¬ 
ported particles. The displacement of a CR particle is expressed by 


ARj — Rj(to): 


( 9 ) 


where the notation j denotes the particle j and R = \/+ y'^ + x = x{t) — x{to). 
The dispersion of displacement R, hereafter a/((A i?)^), is expressed by 


\/«ARp) = 


\ 


^ptls 

E (AR,] 

i=i 


N, 


ptls 


( 10 ) 


where NpUs denotes the total number of particles. 

Here, we show the transport property of particles in our simulations as a function of 
logt, and analyze the diffusion in configuration space. The lines in Figure H] represent 
Y^((AR)2y for protons with energies of Eq = eV, 10^'* eV, and ® eV, corre¬ 

sponding to Rg (gyroradius) > Ax (grid size) (see Figure [2|). The values are averaged 
over Nptis = 10^ particles for each set of energies (eq). We calculated up to t = 10^ 
years. To study the effect of changing the initial location, we defined several boxes for 
the initial CR particle locations. Table [T] shows the physical values of the background 
MHD data in the specific regions, defined as three-dimensional boxes. For the study of 
diffusion in a multiphase medium, we describe the particle distribution and background 
information for the case of a box of size (0.3 pc)^ corresponding to region-I, as given in 
Table [TJ 

The evolution of ((Ai?)^) in Figure [3] is described by standard diffusion, the so- 
called Bohm diffusion, that is the slowest process to be proportional to for all energies 
in our simulations, Eq = eV, 10^"^ eV, and eV. The black dotted line in Figure 

Elis the slope For efficient acceleration of CRs, small diffusion coefficient is required. 
Thus, our results can demonstrate that the diffusion process of CRs in SNRs is can be 
approximated by Bohm diffusion, which may support the assumptions in 271. 


3.1.3. Comparison with Bohm Diffusion Rate 

In this Section, we quantify the diffusion coefficient of CR particles. According to 
the standard assumption of the theory of DSA, accelerated particles diffuse out in the 
plasma, where the diffusive flux obeys Tick’s law: the flux of particles is proportional to 
the gradient of flux density [i^ . It is based on the assumption that charged particles are 
diffuse isotropically by turbulent electromagnetic fields. 

In general, the diffusion coefficient can be expressed as 

D{r,p) = ^Xv, (11) 

where A and v are defined as the mean scattering length and the particle velocity, respec¬ 
tively. The frequently used choice for the diffusion coefficient in DSA is that of Bohm 

6 











diffusion, expressed as 


B 

SB 


2 


( 12 ) 


A 


Ra 


This depends on the degree of turbulence ijil. 

In our MHD simulation, however, the field strength shows the relation SB ~ B, as 
shown in Table 1. Therefore we simply use a simplified expression for Bohm diffusion 
with A ~ i?g: 


R BohmiS'^ P) — 3 ^^^^ 


(13) 


This model, DBohmir, p), has been used by previous work (e.g., [46|'). 

In order to quantify the diffusion of CR particles in our background field, we define 


Rsim (i') 


((AR)^) 

6t 


(14) 


We compare this quantity to the Bohm diffusion coefficient. If ((Ai?)^) is propor¬ 
tional to t, Dsim{t) is constant in time and can be used as the diffusion coefficient in 
the standard diffusion equation for CR particles. Strictly speaking, Dsim{t) is not the 
diffusion coefficient if it is not constant in time. Even in that case, however, Dsim{t) is 
a measure of the diffusion of CR particles. 

Figure shows the ratio of Dgim/DBohm for different values of eo calculated with 
Equations USD and HU)- The values are calculated for 10^ particles at five different 
values of Eq at intervals of logEo = 0-5 eV. Black open and red filled circles correspond 
to < = 10 years and 10^ years, respectively. Both diffusion coefficients, Dgim and DBohm, 
are summarized in Tabled The value of Dsim/DBohm decreases with increasing Eq and 
remains of order unity, 13.5 eV < log Eq < 14.5 eV owing to Bohm diffusion process in 
our simulation. Thus, the Bohm diffusion coefficient description for CR diffusion can be 
justified by our test particle simulations, at least for the shocked SNRs. 


3.1.4- Effect of Electric Field 

As we mentioned in Section 2.2, the electric field in the comoving frame of the ISM 
vanishes in the ideal MHD limit. The strength of the electric field in the computation 
frame is smaller than the magnetic field strength by a factor of w/c (see Equation ([T])). 
Therefore, the force due to the electric field is small in comparison to the force due to 
magnetic field. In this Section, we study the effect of electric field on the diffusion of 
CRs in the ISM. Figure [ 0 ] shows the results of diffusion calculations with and without the 
background plasma {E). It shows that inclusion of E results in slightly larger diffusion 
but the difference is rather limited (< 34%). Therefore, electric fields are less important 
for studying diffusion in configuration space, but they are critical to studying the diffusion 
in energy space. 


3.2. Diffusion Process in a One-phase Medium 

From this Section, we describe the transport properties of particles in configuration 
space for a one-phase medium as a function of log t, and analyze the diffusion in config¬ 
uration space. The values are averaged over Nptis = 10^ particles for each set of energies 
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(eo)- We calculated up to t = 10^ years. We mainly demonstrate the particle distribu¬ 
tion and background information of region-I as given in Table [31 This region is chosen 
to be in the dense, strong magnetic field post-shock region. The calculation method and 
conditions for a one-phase medium are the same as a multiphase medium of Section 3.1. 

3.2.1. Dijfusion in Configuration Space 

The lines in Figure [5] represent a/((A i?)^) for protons with energies of eo = 
eV, 10^^ eV, and ® eV, corresponding to Rg > Ax and Lmin < Rg < Lmax (see 
the gyroradius in Figure 0 and [36[). The uppermost black dashed line in the figure 
represents the slope We observed that the evolution of a/ ((Ai?)^) for Eq = 10^^ eV 
(blue solid line) and eV (red solid line) can be fitted by corresponding to fast 

diffusion, so-called “superdiffusion". Superdiffusion has been studied by S. Xu & H. Yan 
(2013) [dlj and A. Lazarian & H. Yan (2014) [i^. One of the reasons for superdiffusion 
can be the wandering of magnetic field lines in a turbulent medium (see A. Lazarian & H. 
Yan (2014) for details). We demonstrate about superdiffusion property in next Section 
3.2.3. In Figure 8, a/ {{AR)^) for eo = 10^^'^ eV shows standard diffusion, corresponding 
to a simple random walk, represented by black dotted line in Figure [5] This random 
walk can be characterized by a diffusion coefficient, as we will discuss in the Section 3.1. 

3.2.2. Comparison with Rohm Dijfusion Rate 

Figure ini shows the ratio of Dsim/DBohm for different energies (eo), calculated using 
Equations 031) and (see Section 3.1 for review). Applying the Bohm diffusion 
process, the differences in Dsim/DBohm are less than a factor of 5 with energy range 
with 13.5 eV < log eo < 14.5 eV, at both times in Figure [H The differences in the ratio 
of Dsim/DBohm between t = 10 years and 10^ years have grown more than those in Figure 
[3] due to superdiffusive process. However, because the differences in Dsim/DBohm are 
significantly small, thus our result precisely justifies the Bohm diffusion. The diffusion 
coefficients for a one-phase medium are summarized in Table (3] 


3.2.3. Superdiffusion 

Superdiffusion has been reported in simulations by 4^, 411 that is expected to occur 
below the injection scale, (see the injection scale of one-phase medium case in Figure 
d]). We compare our results with S. Xu & H. Yan (2013) to investigate the existence of 
superdiffusion in our simulations. The pre-shock magnetic field lines in our simulations 
are laminar and not turbulent, and the turbulence is created by the RMI driven by 
the propagation of a blast wave. In contrast, S. Xu & H. Yan (2013) drove turbulence 
by solenoidal forcing in Fourier space. 

We calculated the dispersion a/((A^jJ^ as a function of (|A^|| |) up to 10^ years, the 
result of which is shown in Figure [TUI where is perpendicular to the magnetic field 
(B), and ^|| is the displacement along the magnetic field. The upper black dashed lines 
represent a/((A^_l)2 oc respectively. 

The slope close to 1.5 may a manifestation of wandering magnetic field lines analogous 
to Richardson diffusion as discussed in A. Lazarian & H. Yan 2013. This may the reason 
for the superdiffusion found in FigurejSl However, the larger energy particles do not show 
the 1.5 slope of Figure [Till because the gyroradius of larger energy particles is larger than 
the injection scale of turbulence. 
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4. Summary and Discussion 


Young supernova remnants are accompanied by collisionless shock waves in which 
DSA is expected to occur. Gamma rays are emitted by pion decay caused by CR protons 
accelerated in SNRs interacting with molecular clouds. Observations of gamma rays from 
the vicinity of SNRs have shown strong evidence that galactic CR protons are accelerated 


by the shock waves of the SNRs [10|. In DSA theory, efficient acceleration of CRs 
requires a small diffusion coefficient. Indeed, the time-variability seen in synchrotron X- 
ray observations of RX J1713.7-3946 seems to need the slowest diffusion process Q. T. 
Inoue et al. (2^ adopted Bohm diffusion as a diffusion of CRs in the same SNR to explain 
both the gamma-ray spectrum and X-ray variability with magnetic field amplification 
up to mC. In previous research, Bohm diffusion is discussed by considering {B^ /6 
1) (see. Equation (fT^ . 


and [38|). 


In the present study, we performed test particle simulations that describe the propaga¬ 
tion of high-energy CR particles in SNRs in the early evolutionary phase. We investigate 
the diffusion of CR particles in a background medium that is created from numerical sim¬ 
ulation of a realistic interstellar medium. Our results show that the Richtmeyer-Meshkov 
instability can provide enough turbulence downstream of the shock to make the diffu¬ 
sion coefficient close to the Bohm level for energy larger than 30 TeV. Thus, the present 
study indicates that we can use Bohm diffusion coefficient in practical calculations of 
CR particle diffusion, which may support the assumption used in the interpretation of 
high-energy emission from SNRs 


27|. 


It is worth stressing that in this paper we do not deal with the diffusion of particles 
in the medium upstream of the shock. We recall that the DSA mechanism can predict 
the maximum energy of particles needed to explain the observations, only if the diffusion 
coefficient is close to the Bohm one also in the upstream. 

We also identified superdiffusion in a one-phase medium, in which the travel distance 
is proportional to Because the period of superdiffusion is limited in time, the actual 
travel distance is not so different from that described by a Bohm diffusion coefficient, at 
least for particles with energies above 30 TeV. 
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region-1 

post-shock region 

(B) (G) 

2.36e-5 

2.13e-5 

(G) 

3.36e-5 

3.53e-5 

(p) (cm 3) 

2.15 

5.65 

y^iSp)"^ (cm-3) 

1.84 

5.22 


Table 1: Physical values in different simulation regions for the multiphase medium. “region-I” and “post¬ 
shock region" indicate the volume sizes occupied by the initial positions of particles in our simulation. 
region-I and post-shock region have volumes (0.3 pc)^ for the range 1.5 - 1.8 pc in xyz configuration 
space and (2 pc)^ in 1.6 - 1.8 pc of x, 0.7 - 0.9 pc of y, and 1.4 - 1.6 pc of ^ configuration space. In this 
paper, we show the results using region-I for a multiphase medium. 



t = 

10 years 

t = 

100 years 

logeo (eV) 

Dsim (cm^/s) 

Dsohm (cm^/s) 

Dsim (cm^/s) 

Dsohm (cm^/s) 

15.0 

5.839e-H27 

8.000e-h27 

7.520e-h27 

6.584e-h27 

14.5 

2.044e-H27 

2.497e-h27 

2.560e-h27 

2.163e-H27 

14.0 

9.054e-h26 

7.906e-h26 

8.116e-h26 

6.437e-h26 

13.5 

4.745e-h26 

3.193e-h26 

3.193e-h26 

2.038e-h26 

13.0 

3.314e-h26 

2.015e-h26 

2.015e-h26 

6.225e-h25 


Table 2: Differences between Dsim vs. a function of £o- We used the region-I volume results 

shown in Tabled 


region-I region-II post-shock region 


(B) (G) 1.262e-5 6.033e-6 1.102e-5 

(G) 1.092e-5 9.306e-6 1.123e-5 

(p) (cm-3) 1.303 0.12 1.025 

v/(Vp)2 (cm-3) 1.017 0.313 0.917 


Table 3: Physical values in different simulation regions for one-phase medium. “region-I”, “region-II”, 
and “post-shock region" indicate the volume sizes occupied by the initial positions of particles in our 
simulation. region-I and region-II have volumes of (0.3 pc)^ for the range 0.7 - 1 pc in xyz configuration 
space of our simulated data and (0.2 pc)^ for 0.8 - 1.0 pc in x, 0.7 - 0.9 pc in y, and 1.0 - 1.2 pc in 2 
configuration space of our simulated data. These regions are chosen to be in the dense, strong magnetic 
field post-shock region. Our simulation is performed in region-I for a one-phase medium. 



t = 

10 years 

t = 

100 years 

logeo (eV) 

Dsim (cm^/s) 

Dsohm (cm^/s) 

Dsim (cm^/s) 

Dsohm (cm^/s) 

15.0 

1.120e-h28 

1.334e-h28 

2.575e-H28 

1.277e-H28 

14.5 

3.593e-h27 

3.073e-h27 

1.189e-H28 

4.134e-H27 

14.0 

1.318e-h27 

7.920e-h26 

2.916e-h27 

9.189e-h26 

13.5 

7.578e-h26 

2.126e-h26 

8.687e-h26 

2.071e-h26 

13.0 

6.374e-h26 

6.191e-h25 

7.136e-h26 

4.740e-h25 


Table 4: Difference between Dsim and as a function of so¬ 


il 




















T, Inoue et al. (2012) 



From the data of T. Inoue et al. (2013) 
one-phase medium 



Figure 1: Fourier power spectrum of SNR turbulence. Figure on left and right indicate for multiphase 
and one-phase medium, respectively. The upper thin solid line in figure indicates the velocity field 
and dashed line indicates the magnetic field. The uppermost thick solid line in figure represents the 
Kolmogorov law 



13 



12 

slope=0.5 

> 


/ 


11 


Cj? 


(' 


10 

u 



-1.5 -1 -0.5 0 0.5 1 1.5 2 



log t (year) 


log t (year) 


Figure 2: The magnitude of the change in particle energy, Ae, as a function of logt for protons with an 
initial energy of eg = 10^ TeV (left) and eg = 10^ TeV (right). We show several cases (solid line) for 
this Figure The uppermost dotted line in figure shows the slope 
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Figure 3: The gyroradius (Rg) in units of pc as a function of cosmic ray initial energy (so)- The results 
at 10 years and 10^ years are shown by black open and filled circles, respectively. Red star symbol is for 
the case t = 0. 


13 







- 0.5 0 0.5 1 1.5 2 

log t (year) 


Figure 4: Travel distance ( yj {(Ai?)^) ) of cosmic ray particles (so) as a function of log t. Green, blue, 
and red lines indicate the curves of eo = eV, 10^“^ eV and eV, respectively. All lines are 

similar to the uppermost black dashed line that is proportional to t^'^, corresponding to the slowest 
process. The travel distances with -0.5 < log i < 0.5 are proportional to t. 
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Figure 5: The diffusion coefficient as a function of cosmic ray energy (eo). Vertical axis shows the ratio 
between the simulation values {Dsim) the estimated values {DBohm) from Equations ns and d, 
shown as a function of log so by solid lines. The results at 10 years and 10^ years are shown by open 
black and filled red circles, respectively. 
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Figure 6: Comparison between ratios of results with E and without E for the cosmic ray energy range 
(eo) 10 TeV to 10^ TeV. 
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Figure 7: The gyroradius (Rg) in units of pc as a function of cosmic ray initial energy (so)- The results 
at 10 years and 10^ years are shown by black open circles and filled circles, respectively. Red star symbol 
is for the case t = 0. 



17 







-0.5 0 0.5 1 1.5 2 

log t (year) 


Figure 8: Travel distance ( ((Ai?)^) ) of cosmic ray particles (eo) as a function of log t. Green, blue, 

and red lines indicate the curves of so = eV, 10^"^ eV and eV, respectively. The red and 

blue lines are similar to the uppermost black dashed line that is proportional to corresponding to 

the fast process (superdiffusive). And green line is similar to the uppermost black dashed line that is 
proportional to t^-^, corresponding to the slowest process. 
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Figure 9: The gyroradius and the diffusion coefficient as a function of cosmic ray energy. The left panel 
is gyroradius in units of pc, shown by dashed lines. The right panel is the ratio between the simulation 
values (Dsim) ^nd the estimated values (-D^ohm) from Equations i(T3) i and (11411 . shown as a function 
of log£o by solid lines. The results at 10 years and 10^ years are shown by open black and filled black 
circles, respectively. Stars are for the case t=0. 
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Figure 10: Relation between travel distances in the directions parallel and perpendicular to the magnetic 
field lines. Green, blue, and red lines indicate the curves of eg = eV, 10^^ eV and eV, 

respectively. The black dotted line describes oc (|A^|||)^’^ in limited time. 
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